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Abstract 

A  two-phase  plume  flow  from  a  small  aluminized  pro¬ 
pellant  side  thruster  interacting  with  rarefied  atmo¬ 
sphere  at  120  km  has  been  examined  numerically. 
A  three  step  continuum-kinetic  approach  has  been 
used,  with  the  Navier-Stokes  equations  solved  inside 
the  nozzle,  and  a  2D/3D  DSMC  method  employed 
to  compute  the  plume  nearfield  and  then  the  plume- 
atmosphere  interaction  region.  At  each  of  these  steps, 
a  two-way  gas-particulate  coupling  has  been  used. 
The  DSMC  implementation  uses  molecular  fluxes  to 
calculate  the  number  of  gas-particulate  collisions,  and 
is  based  on  a  statistical  approach  to  calculate  deflec¬ 
tion  angles.  A  sensitivity  study  of  various  parameters 
of  the  approach  is  performed.  The  importance  of  two- 
way  coupling,  particle  radiative  cooling,  and  molecule 
accommodation  on  particle  surface  are  analyzed. 

1  Introduction 

The  combustion  and  relaxation  processes  in  rocket 
propulsion  systems  result  in  the  formation  of 
particulates  of  different  size  and  nature.  Alu¬ 
minum  oxide  particles,  associated  with  all  types 
of  aluminized  propellants,  compose  up  to  30%  of 
the  exhaust  mass  fraction  and  may  significantly 
impact  both  the  gas  flow  inside  the  nozzle  and 
plume- atmosphere  interaction  phenomena.  Unburnt 
drops  of  liquid  propellants  are  another  important 
class  of  particles.  Although  large  in  size,  they  do 
not  affect  the  bulk  gas  flow  due  to  their  small 
mass  fractions;  they  may,  however,  cause  signifi¬ 
cant  contamination  problems.  Soot  must  also  be 
considered  for  liquid  propellant  thrusters;  these  par¬ 
ticulates  are  not  expected  to  significantly  impact  the 
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gas  flow  properties,  but  are  important  contributors  to 
plume  radiation  signature. 

Significant  efforts  have  been  made  in  the  aerospace 
community  toward  accurate  prediction  of  aluminum 
oxide  particle  impact  on  thruster  performance  and  ex¬ 
haust  plume  structure  (see,  for  example,  [1]).  Many 
of  the  numerical  efforts  however  were  limited  to  low- 
altitude,  high  ambient  density  flow  regimes.  The 
continuum  CFD  based  approaches  can  not  be  used 
to  model  the  influence  of  particulates  on  the  plume- 
atmosphere  interaction  at  high  altitudes,  where  the 
flow  nonequilibrium  is  too  strong  and  these  ap¬ 
proaches  become  unsuitable.  A  kinetic  approach, 
such  as  the  direct  simulation  Monte  Carlo  (DSMC) 
method,  has  to  be  used  to  obtain  credible  informa¬ 
tion  on  these  flows. 

The  main  objectives  of  this  work  are  to  com¬ 
pare  gas-only  and  two-phase  steady-state  three- 
dimensional  DSMC  predictions  of  the  flowfield  sur¬ 
rounding  a  small  side-mounted  solid  propellant  at¬ 
titude  control  thruster  at  high  altitudes,  and  to  as¬ 
sess  the  impact  of  particles  on  shock  layer  structure 
and  plume  observables.  To  the  best  of  our  knowl¬ 
edge,  this  work  is  the  first  application  of  the  DSMC 
method  to  model  full  3D  two-phase  plume-atmosphere 
interaction.  Special  attention  is  paid  to  the  influence 
of  different  simulation  parameters  and  flow  phenom¬ 
ena,  such  as  one-  or  two-way  coupling,  gas-particulate 
interactions,  and  particle  radiation  cooling,  on  gas 
and  particulate  density  and  temperature  distribu¬ 
tions.  The  present  study  focuses  on  a  flow  from  a 
130  N  aluminized  propellant  thruster  into  atmosphere 
at  an  altitude  of  120  and  a  free  stream  speed  of  5  km/s. 

Large  particle  loadings  for  the  cases  under  consider¬ 
ation  prevent  the  application  of  a  simple  overlay  par¬ 
ticle  tracking  approach  and  a  two-way  coupling  needs 
to  be  used.  The  first  implementation  of  a  two-way 
coupling  in  the  DSMC  method  was  presented  in  [2]. 
The  approach  proposed  in  this  work  is  different  from 
[2],  and  the  details  of  the  present  approach  are  given 
below.  Similar  to  [2],  the  model  [3]  was  used  to  sim- 
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ulate  the  impact  of  gas  molecules  on  aluminum  oxide 
particulates.  The  current  algorithm  does  not  use  an 
analytic  approximation  for  the  molecule  deflection  an¬ 
gle,  but  implements  an  approach  based  on  a  statistical 
determination  of  the  reflection  point,  and  reflection 
procedure  in  accordance  with  the  Maxwell  model  with 
specified  tangential  and  energy  accommodation  coeffi¬ 
cients.  The  molecule-particle  collision  rate  is  based  on 
the  particle  number  flux  values,  and  not  on  the  NTC 
scheme.  All  energies  (translational  and  internal)  of 
molecules  are  included  in  the  energy  transfer,  and  the 
Maxwell  diffuse-specular  model  is  used  to  calculate 
after-collision  molecular  states.  The  proposed  proce¬ 
dure  has  linear  dependence  of  computational  time  on 
the  number  of  molecules  and  particles. 

2  General  modeling  strategy 

The  modeling  of  the  plume-atmosphere  interaction  is 
performed  in  three  steps.  First,  the  VIPER  two-phase 
CFD  code  [4,  5]  is  utilized  to  obtain  the  solutions  in¬ 
side  the  nozzle.  The  VIPER  solution  at  the  nozzle  exit 
is  then  used  as  the  inflow  boundary  condition  for  an 
axisymmetric  DSMC  computation  of  the  plume  near 
field.  As  a  result  of  the  second  step,  a  starting  surface 
for  a  subsequent  3D  DSMC  computation  is  available 
at  a  distance  from  the  exit  were  the  plume  density  is 
too  high  for  the  ambient  atmosphere  to  have  any  no¬ 
ticeable  impact  on  the  plume.  The  third  step  is  full 
3D  DSMC  modeling  of  the  plume  -  atmosphere  inter¬ 
action  that  also  includes  chemical  reactions  between 
gas  species.  SMILE  [6]  DSMC-based  computational 
tool  is  used  for  the  second  and  third  steps,  extended 
to  include  particulates  coupled  with  the  gas  flow.  All 
three  stages  incorporate  a  two-way  coupling  between 
gas  molecules  and  alumina  particles. 

2.1  VIPER  code 

VIPER[4,  5]  is  an  axisymmetric  Parabolized  Navier- 
Stokes  (PNS)  code  that  includes  finite  rate  gas  chem¬ 
istry,  multiphase  capability  (via  a  two-way  coupled 
Lagrangian  method),  and  a  variety  of  mostly  empiri¬ 
cal  models  for  gas-particulate  interaction  and  partic¬ 
ulate  evolution  phenomena.  The  PNS  scheme  is  ap¬ 
plied  from  the  sonic  line  near  the  throat  to  the  nozzle 
exit  plane.  Separate  methods  are  provided  to  model 
the  combustion  chamber  and  converging  section.  The 
combustion  chamber  pressure  and  temperature  were 
assumed  to  be  3.85  atm  and  2900  K,  respectively. 
These  conditions  result  in  an  exit  plane  gas  pressure 
of  approximately  one  percent  of  one  atmosphere,  and 
are  essentially  the  lower  limits  for  which  a  successful 
VIPER  run  could  be  made.  Under  these  conditions, 
it  is  expected  that  the  VIPER  results  give  a  quali- 
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tative,  rather  than  quantitative,  prediction  of  the  2- 
phase  nozzle  flow,  suitable  for  engineering  studies. 

2.2  SMILE  code 

The  2D/3D  DSMC-based  code  SMILE  [6]  has  been 
used  as  the  principal  kinetic  approach.  The  impor¬ 
tant  features  of  SMILE  that  are  relevant  to  this  work 
are  parallel  capability,  different  collision  and  macropa¬ 
rameter  grids  with  manual  and  automatic  adapta¬ 
tions,  and  spatial  weighting  for  axisymmetric  flows. 

The  major  ant  frequency  scheme  [7]  was  used  to  cal¬ 
culate  inter  molecular  interactions.  The  intermolecu- 
lar  potential  was  assumed  to  be  a  variable  hard  sphere 
[8] .  Energy  redistribution  between  the  rotational  and 
translational  modes  was  performed  in  accordance  with 
the  Larsen-Borgnakke  model.  The  total  collision  en¬ 
ergy  model  [8]  was  used  to  model  chemical  reactions  in 
gas.  Species  weights  were  used  for  particulate  species. 
A  temperature-dependent  rotational  relaxation  num¬ 
ber  was  used.  The  reflection  of  molecules  on  the  noz¬ 
zle  and  rocket  surface  was  assumed  to  be  diffuse  with 
complete  energy  and  momentum  accommodation.  An 
extension  of  SMILE  to  simulate  two-phase  flows  is  dis¬ 
cussed  in  the  following  section. 

3  DSMC  implementation  de¬ 
tails 

The  present  DSMC  implementation  uses  the  following 
assumptions:  (i)  particulates  are  spherical,  and  their 
rotation  does  not  impact  the  flow,  (ii)  particulates  are 
uniform,  and  have  no  internal  temperature  gradients, 
and  (iii)  the  gas  mean  free  path  is  larger  than  the 
particle  diameter,  so  the  gas  flow  is  essentially  free- 
molecular  with  regard  to  particle  sizes.  This  means 
that  there  are  no  significant  gas  gradients  near  the 
particle  surface. 

The  most  important  effects  considered  in  the  model 
are  particle  drag  and  radiative  and  collisional  heat¬ 
ing/cooling  of  particles,  and  impact  of  particles  on 
gas  (two-way  coupling).  The  following  effects  are  not 
considered  at  this  time:  chemical  reactions  on  the  sur¬ 
face,  processes  that  occur  in  liquid  droplets  (evapo¬ 
ration/condensation,  coalescence),  and  particle  phase 
change. 

The  most  important  features  of  the  current  imple¬ 
mentation  that  distinguish  it  from  [2]  are  as  follows. 
The  implementation  (i)  does  not  account  for  particle 
rotation:  since  particles  in  rocket  plumes  are  largely 
spherical  [9],  particle  rotation  is  considered  unimpor¬ 
tant,  (ii)  uses  a  statistical  approach  to  calculate  de¬ 
flection  angles  of  molecules  on  particle  surface,  (iii) 
gives  flexibility  of  using  arbitrary  model  of  reflection  of 
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molecules  on  particulates,  and  includes  vibrational  en¬ 
ergy  transfer  in  molecule-particle  collisions,  (iv)  uses 
molecular  fluxes  and  not  the  NTC  scheme  to  calculate 
the  number  of  molecule-particle  collisions,  (v)  utilizes 
arbitrary  inflow  distribution  of  particles  in  terms  of 
sizes  and  velocities,  and  (vi)  employs  different  averag¬ 
ing  strategy  for  particle  force  /  energy  transfer  calcu¬ 
lations. 

3.1  Force  and  heat  flux  computation 

The  force  F  and  heat  flux  Q  on  each  particle  exerted 
by  gas  molecules  are  calculated,  similar  to  [3],  as  the 
sum  over  individual  force  and  heat  flux  contributions 
of  all  surrounding  molecules  in  a  given  cell, 

F  =  J  Fs(c)f(c)dc,  Q  =  J  Qs{c)f{c)dc, 

where  c  =  u  —  up,  u  is  the  molecular  velocity  vector, 
and  Up  is  the  particle  velocity  vector.  This  is  approxi¬ 
mated  by  the  summation  over  all  molecules  i  in  a  cell 
at  a  given  timestep, 

N  N 

i= 1  i= 1 


degrees  of  freedom,  v  denotes  summation  over  vibra¬ 
tional  modes,  and  d  is  the  mode  degeneracy. 

Force  and  heat  transfer  computations  use  the  aver¬ 
age  properties  u  and  Tp  for  particulates  that  depend 
on  particulate  species,  diameter,  and  spatial  cell.  Ev¬ 
ery  time  step  the  particle  velocities  are  recalculated 
as  follows, 


v(t  +  At)  =  v(t)  +  F7TRpAt/mp , 

where  v  and  mp  are  particulate  velocity  vector  and 
mass.  F  is  the  force  per  unit  area  (see  above)  multi¬ 
plied  by  the  actual  area  of  the  particulate  irRp.  Par¬ 
ticle  temperature  is  recalculated  as 


T(t  +  At)  =  T(t)  +  AtirRp/ (mpCp)Q+ 


Pp  Rp  Cj 


-HT*as 


Here,  Cp  is  the  specific  heat  capacity  of  particulates,  Q 
is  the  heat  flux  per  unit  area  calculated  at  the  previous 
step,  a  is  the  Stefan-Boltzmann  constant,  and  Tgas 
is  the  temperature  of  the  surrounding  gas  (the  cell 
temperature) .  The  second  term  in  the  right  hand  side 
of  the  equation  corresponds  to  the  radiation  gain  or 
loss.  Note  that  an  emissivity  of  one  is  assumed  here. 


It  is  also  possible  to  compute  the  average  force  and 
heat  flux  over  M  consecutive  time  steps  as, 


F 


M  N 


j= 1 i= 1 


Introducing  the  momentum  accommodation  coeffi¬ 
cient  am  and  the  translational,  rotational,  and  vibra¬ 
tional  energy  accommodation  coefficients  at ,  ar,  and 
av ,  one  can  write  the  expression  for  the  force  per  unit 
area, 


Fi 


=  m 


-  C- 

Vcell 


(1  +  gam(l!  —  OLt))g  + 


2 -Tp 
m 


3.2  Two-way  coupling 

To  account  for  the  impact  of  particulates  on  gas,  one 
needs  to  model  the  particle-molecule  collision  pro¬ 
cess.  Two  principal  issues  need  to  be  addressed  in 
this  respect,  (i)  computation  of  the  number  of  gas- 
particulate  collisions,  and  (ii)  modeling  of  each  indi¬ 
vidual  collision  (change  molecular  velocities  and  ener¬ 
gies). 

The  total  number  of  collisions  of  gas  molecules  with 
a  specific  particle  is  calculated  as  the  free-molecular 
number  flux  of  molecules  to  the  particle  surface,  which 
may  be  written  as 


where  m  is  the  molecular  mass,  Fnum  is  the  ratio  of 
real  to  simulated  molecules,  Vceu  is  the  cell  volume, 
Rp  is  the  particle  radius,  g  =  |c|,  and  Tp  is  the  particle 
temperature. 

For  the  heat  flux  per  unit  area, 

Qi  =  am  ~  2 kTp)+ 


ar(Erot  —  ^-—kTp)  +av(EVib  — 


d  k  0„ 


exp  {9V/Tp) 


7>}- 


Here,  Erot  and  are  the  molecular  rotational  and 
vibrational  energies,  f rot  is  the  number  of  rotational 


where  n  is  the  gas  number  density,  c  is  the  molecular 
velocity  with  respect  to  particle  velocity,  and  D  is 
particle  diameter. 

In  this  work,  the  majorant  collision  frequency  [7] 
methodology  has  been  applied,  so  that  the  number  of 
majorant  collisions  is 


Nr.. 


,nc 

=  <T’D > 


max-) 


and  the  acceptance-rejection  technique  is  used  to  se¬ 
lect  physical  collisions. 

Let  us  now  describe  the  algorithm  for  molecule- 
particle  collisions.  It  has  two  parts,  (a)  obtain  the 
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point  where  molecule  hits  particle  surface,  and  (b) 
perform  reflection.  In  part  (a)  the  surface  reflection 
point  is  selecteded  stochastically.  After  the  reflection 
point  is  determined,  the  normal  to  the  surface  may 
found  assuming  that  particles  are  spherical.  Then, 
(b),  new  molecular  states  (velocities  and  internal  en¬ 
ergies)  are  calculated  in  the  same  way  as  for  usual  re¬ 
flection  of  molecules  on  surfaces,  and  using  specified 
accommodation  coefficients.  Note  that  molecular  co¬ 
ordinates  are  not  changed.  The  particle  properties  are 
not  changed  either  since  the  impact  of  gas  on  parti¬ 
cles  is  included  through  the  momentum  and  heat  flux 
modeling  (see  previous  sections). 

Part  (b)  is  straightforward  and  therefore  will  not 
be  described  here.  Part  (a)  consists  of  the  following 
steps. 

1.  Calculate  the  reflection  point  in  a  local  coordi¬ 
nate  system.  The  local  coordinate  system  is  set 
so  that  the  relative  velocity  vector  is  parallel  to 
the  X  axis.  The  system  is  moving  with  the  parti¬ 
cle,  so  that  the  particle  is  still,  and  the  molecule 
approaches  with  a  negative  velocity  (see  Fig.  1). 
The  following  approximation  is  used  here, 

R  Rmol  T  Rpart  Rpart' 

The  impact  parameter  b  is  selected  randomly  as 
b  =  RfR. 

The  reflection  point  is  initially  located  on  XY 
plane,  with  its  coordinates  determined  by  a  vec¬ 
tor 


The  equation  of  transformation  is  written  as 


ell 

612  ei3  \  / 

y/R2  -  b2 

^21 

622  623  1*1 

bcos<p 

631 

632  633  /  \ 

b  sirup 

where 

6n  =  vX7  ei2  =  —vy,  ei3  =  —vz, 

vz  ,  vVvz 

e2i  =  vy,  e22  =  — - b  vx,  e23  =  - , 

l  +  vx  l  +  vx 

Vy  V  z  Vy 

e3i  =  vzi  e32  =  - ,  e33  =  - — - +vx. 

l  +  vx  1  +  vx 

4.  Finally,  the  needed  normal  vector  is  determined 
by  the  following  directional  cosines, 

cos  ipi  =  e\\ y/ R2  —  b2  +  e^b  cos  <p  +  6136  sin  0, 

cos  ip 2  =  e2i  y/  R2  —  b2  +  e22^>  cos  <p  +  e23^  sin  <p, 
cos  -03  =  e3i  y/  R2  —  b2  +  e^b  cos  (p  +  e^b  sin  p. 

4  Model  Verification 

Prior  to  being  used  to  model  complex  3D  flows,  the 
extended  SMILE  code  has  been  extensively  tested  to 
verify  its  robustness  and  accuracy.  An  example  of  the 
verification  is  shown  in  Figs.  2  and  3  where  the  veloc¬ 
ity  and  temperature  of  particulates  introduced  to  a 
uniform  gas  stream  are  presented.  The  initial  partic¬ 
ulate  temperature  is  1000  K  and  velocity  is  1000  m/s, 
and  the  corresponding  gas  properties  are  2000  K  and 
2000  m/s.  The  gas  and  particle  number  densities  are 
1022  molecule/m3  and  1012molecule/m3,  respectively. 
Note  that  the  obtained  DSMC  solution  is  in  excellent 
agreement  with  the  analytic  solution. 


2.  Rotate  the  reflection  point  around  X  axis  at  an 
angle  <p  =  2nRf ,  which  will  give  the  coordinate 
of  the  reflection  point  in  the  local  coordinate  sys¬ 
tem, 


£0  =  Axq 


10  0 
0  cos  <p  —  sin  <p 
0  sin  <p  cos  <p 


3.  Transform  the  local  (xi)  to  the  global  coordinate 
system  via  the  directional  cosines, 


_  RQx  _  Rfly  _  Rgz 

COS  Qq  :  — — — ,  COS  0^2  =  -7 — p,  COS  0^3  :  — — — . 

\g\  \g\  \g\ 


5  Flow  Conditions 

The  test  case  investigated  in  this  work  is  a  two- 
phase  flow  from  a  130  N  RCS  aluminized  propellant 
thruster  into  atmosphere  at  an  altitude  of  120  km. 
The  thruster  was  located  at  a  side  of  a  small  rocket 
that  travels  at  a  speed  of  5  km/s  and  zero  incidence. 
A  conical  nozzle  was  used  with  a  throat  diameter 

1.6  cm,  an  exit  diameter  of  8.8  cm,  and  a  half- angle 
of  15  deg.  The  alumina  particle  mass  loading  in  the 
stagnation  chamber  was  assumed  to  be  18%.  The  alu¬ 
mina  particles  were  assumed  to  have  a  diameter  of 

3.6  /im  that  did  not  change  in  the  simulations.  The 
gas  composition  (mass  fraction)  at  the  nozzle  throat 
was  X[CO]=0.5186,  X[HC1]=0.2559,  X[N2]=0.1170, 
X[H2]=0.0479,  X[H20]=0.0468,  X[C02]=0.0138.  The 
nozzle  surface  temperature  was  assumed  to  be  1000  K. 
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The  rocket  was  modeled  a  blunted  cone  cylin¬ 
der,  and  the  thruster  was  positioned  on  the  cylin¬ 
der  immediately  after  the  cylinder-cone  junction.  The 
cone  length  was  1.73  m,  and  the  cylinder  length  and 
thickness  were  2.6  m  and  0.4  m,  respectively.  The 
free  stream  temperature  and  number  density  were 
354  K  and  4.733  molecule/m3,  respectively.  The 
free  stream  composition  was  X[N2]=0.73,  X[O]=0.18, 
X[02]=0.09.  The  reaction  set  used  in  this  work  was 
taken  from  [10]. 

6  Results  and  discussion 

6.1  Flow  Inside  the  Nozzle 

The  computations  of  the  diverging  part  of  the  nozzle 
were  performed  using  both  VIPER  and  SMILE  codes. 
The  throat  conditions  obtained  from  the  VIPER  so¬ 
lution  of  the  converging  part  of  the  nozzle  were  used 
in  the  DSMC  modeling.  Although  over  10  million 
molecules  and  3  million  cells  were  used  in  the  DSMC 
simulations,  these  numbers  are  about  two  orders  of 
magnitude  smaller  than  those  required  to  strictly  sat¬ 
isfy  the  corresponding  requirements  of  the  DSMC 
method  (the  throat-based  Knudsen  number  is  about 
5  •  10-5),  and  the  DSMC  results  should  therefore  be 
considered  as  qualitative. 

Comparison  of  the  continuum  and  kinetic  solutions 
is  presented  in  Fig.  4,  where  the  temperature  fields  are 
shown.  Note  that  the  DSMC  computational  domain 
was  extended  beyond  the  nozzle  exit  ( X  =  0.1  m).  It 
is  seen  that  there  is  a  reasonable  agreement  between 
the  two  solutions,  with  the  larger  difference  observed 
further  from  the  nozzle  centerline.  The  two  solutions 
are  close  near  the  centerline.  It  is  interesting  to  note 
that  this  is  the  region  where  all  the  particulates  are 
concentrated.  Relatively  large  mass  of  particulates 
results  in  their  small  divergence  from  the  nozzle  axis, 
as  illustrated  in  Fig.  5.  The  blue  colored  region  corre¬ 
spond  to  the  part  of  the  flow  where  no  particulates  was 
observed.  The  difference  between  the  continuum  and 
kinetic  solutions  is  primarily  attributed  to  different 
approaches  to  account  for  the  drag  force  and  heat  flux 
on  particles  exerted  by  gas  molecules  (a  continuum 
technique  is  used  in  the  VIPER  code,  whereas  the 
above  free-molecule  approach  is  utilized  in  SMILE). 

The  interaction  between  gas  and  particulates  causes 
a  drop  in  temperature  of  particulates,  from  2,700  K 
at  the  throat  to  about  2,500  K  at  the  nozzle  exit.  The 
change  in  particulate  velocity  is  more  pronounced, 
from  about  700  m/s  at  the  throat  to  about  1,200  m/s 
at  the  exit  (see  Fig.  6). 

Another  important  fact  established  in  the  nozzle 
flow  modeling,  in  addition  to  the  concentration  of  par¬ 
ticulates  in  the  coreflow,  is  their  visible  deviation  from 
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the  nozzle  axis.  This  behavior  is  illustrated  in  Fig.  7, 
which  shows  that  the  density  of  particulates  reaches 
its  maximum  values  at  some  distance  from  the  axis. 
This  is  attributed  to  the  particulate  transport  inside 
the  converging  part  of  the  nozzle,  where  the  nonuni¬ 
formity  in  particle  concentration  is  created.  Particles 
lag  behind  the  gas  flow  in  the  converging  part,  and  as 
the  result  move  closer  to  the  surface,  and  at  the  throat 
plane  form  a  maximum  located  between  the  centerline 
and  the  surface. 

At  high  particle  loading  values  typical  for  alu¬ 
minized  propellant  thruster,  the  gas  influences  the 
particulates,  but  the  latter  may  also  significantly  im¬ 
pact  the  gas  properties.  In  order  to  examine  the  in¬ 
teraction  between  gas  and  particulates,  in  addition  to 
the  baseline  two-way  coupling  modeling  the  DSMC 
computations  were  also  performed  with  a  one-way 
coupling  (particulates  assumed  not  to  affect  the  gas 
molecules) .  Comparison  of  the  gas  axial  velocity  fields 
for  simulations  with  one-way  and  two-way  coupling 
are  presented  in  Fig.  8.  Note  a  qualitative  difference 
between  these  two  cases.  When  the  one-way  coupling 
is  included,  the  gas  velocity  inside  the  nozzle  mono- 
tonically  decreases  with  the  radial  distance  from  the 
nozzle  centerline.  The  interaction  with  particulates 
that  travel  at  lower  speeds,  though,  results  in  a  sig¬ 
nificant  minimum  in  gas  velocity.  The  velocity  then 
increases  at  the  nozzle  axis  due  to  the  local  minimum 
in  particulate  density  in  that  region  (see  Fig.  7). 

Lower  gas  velocities  in  the  coreflow,  in  turn,  result 
in  lower  particulate  velocities  for  the  two-way  coupling 
case  as  compared  to  the  one-way  coupling,  as  shown 
in  Fig.  9.  The  difference  in  the  particulate  velocities 
is  significantly  smaller,  however,  than  for  the  gas  ve¬ 
locities.  At  the  nozzle  throat  it  amounts  to  about 
100  m/s,  versus  over  500  m/s  for  gas  velocities. 

One  of  the  important  parameters  of  the  two-phase 
flow  modeling  is  the  particle  size  and  mass.  Although 
these  parameters  are  relatively  well  known  for  alu¬ 
minized  propellant  thrusters  as  compared  to  liquid 
propellant  engines,  there  is  still  some  uncertainty  in 
particulate  properties.  To  study  the  sensitivity  of  flow 
results  to  the  particulate  mass,  a  DSMC  computation 
has  also  been  performed  for  pure  aluminum  particles, 
whose  mass  is  about  1.5  times  smaller  than  for  the 
baseline  alumina  case. 

The  general  impact  of  the  decrease  in  the  mass  of 
particles  is  their  visibly  larger  divergence  from  the 
centerline  due  to  a  larger  influence  of  the  drag  force. 
This  effect  is  shown  in  Fig.  10  where  the  particle  ax¬ 
ial  velocity  fields  are  given  for  the  two  cases.  Note 
also  that  aluminum  particles  accelerate  to  speeds  up 
to  1,300  m/s  at  the  nozzle  exit  versus  about  1,200  m/s 
for  alumina  particles. 

There  is  also  a  noticeable  impact  of  the  particle 
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mass  on  gas  properties,  illustrated  in  Fig.  11  for  the 
axial  velocity.  As  expected,  the  gas  velocity  in  the 
region  of  particle-gas  interaction  is  lower  for  the  heav¬ 
ier  particles.  The  particle-free  region  is  not  affected, 
though. 

The  gas-particle  interaction  may  significantly  im¬ 
pact  the  infrared  radiation  signatures,  and  one  of  the 
primary  modeling  issues  is  the  accurate  prediction 
of  vibrational  temperatures  of  different  species.  The 
main  mechanism  of  influence  of  gas-particle  collisions 
on  vibrational  temperature  is  the  accommodation  of 
vibrational  energy  on  particulate  surface.  To  examine 
the  possible  impact  of  vibrational  energy  accommo¬ 
dation,  computations  where  performed  for  the  com¬ 
plete  (baseline)  and  no  vibrational  energy  accommo¬ 
dation.  The  results  for  the  vibrational  temperature 
of  one  of  the  gas  species  are  presented  in  Fig.  12. 
The  computations  show  that  the  vibrational  energy 
accommodation  effect  on  vibrational  temperature  is 
relatively  small,  amounting  to  about  10%  at  the  noz¬ 
zle  exit  plane. 

6.2  Modeling  of  the  Plume  Near  Field 

The  macroparameters  at  the  nozzle  exit  plane  ob¬ 
tained  with  the  VIPER  code  were  used  as  the  bound¬ 
ary  condition  for  the  plume  near  field  calculations. 
The  kinetic  approach  was  employed  for  this  stage,  and 
about  4  million  molecules  and  2  million  cells  used  in 
the  computations  have  provided  the  required  accuracy 
and  resolution  for  the  DSMC  method.  The  computa¬ 
tional  domain  extended  2.5  m  downstream  from  the 
nozzle  exit,  where  the  plume  is  still  dense  enough  to 
neglect  the  impact  of  the  free  stream  on  the  coreflow. 
Only  plume  species  have  therefore  been  considered  in 
the  simulation,  with  no  chemical  reactions  included. 

Although  this  stage  has  been  used  primarily  to  pro¬ 
vide  a  starting  surface  for  the  subsequent  3D  compu¬ 
tations,  it  has  also  allowed  us  to  examine  the  impact  of 
the  two-way  coupling  that  has  been  previously  shown 
important  for  the  flow  inside  the  nozzle.  The  gas  ax¬ 
ial  velocity  profiles  are  given  in  Fig  13  for  the  one-way 
and  two-way  coupling  cases.  Note  that  the  nozzle  exit 
plane  is  located  at  X=0.1  m.  As  expected,  there  is  no 
visible  impact  of  particulates  on  gas  in  the  particle- 
free  regions  of  the  plume.  In  the  coreflow,  where  the 
concentration  of  particles  is  considerable,  the  two-way 
coupling  results  in  a  decrease  of  gas  velocities  by  a  few 
percent. 

Although  the  influence  of  gas-particle  collisions  on 
gas  is  smaller  in  the  plume  than  inside  the  nozzle,  this 
effect  is  still  noticeable,  and  is  most  pronounced  for 
the  gas  translational  temperature  fields  (see  Fig.  14). 
The  difference  originates  in  the  first  few  centimeters 
downstream  from  the  nozzle  exit,  where  the  concen¬ 
tration  of  particulates  is  maximum. 
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The  influence  of  two-way  coupling  on  particulate 
properties  is  very  small,  and  the  particle  flow  fields 
obtained  by  one-way  and  two-way  coupling  are  prac¬ 
tically  identical.  The  particle  temperature  field  in 
the  plume  near  field  is  shown  in  Fig.  15.  It  is  seen 
that  the  particle  temperature  decreases  downstream; 
this  decrease  is  attributed  to  the  radiative  cooling  and 
not  the  gas-particle  collisions.  The  computation  with 
the  radiative  cooling  turned  off  has  shown  no  visible 
change  in  particle  temperature.  The  computations 
also  showed  that  the  particle  velocities  do  not  change 
along  the  stream  lines,  proving  therefore  that  the  in¬ 
fluence  of  the  particle  drag  in  the  nearfield  is  negligi¬ 
ble. 

6.3  Plume- Atmosphere  Interaction 

The  3D  modeling  of  the  plume-atmosphere  interaction 
was  performed  with  the  SMILE  code  using  a  starting 
surface  obtained  using  the  plume  near  field  solution 
described  in  the  previous  section.  The  starting  surface 
was  generated  along  the  density  isoline  with  a  value  of 
2  •  1021  kg/m3;  it  expands  about  0.45  m  downstream 
from  the  nozzle  exit.  An  elliptic  distribution  function 
was  used  for  the  plume  molecules  that  enter  the  com¬ 
putational  domain  from  that  starting  surface.  The 
size  of  the  computational  domain  is  10.4  m  in  the  X 
direction,  7.4  m  in  the  Y  direction,  and  2.6  m  in  the  Z 
direction  (a  symmetry  plane  was  used  at  Z=0).  Here, 
X  coincide  with  the  plume  axis,  and  Y  is  parallel  to 
the  free  stream  vector.  A  total  of  about  10  million 
simulated  molecules  and  5  million  collision  cells  was 
used  in  the  3D  computations. 

The  sensitivity  studies  have  shown  that  the  impact 
of  gas-particle  collisions  is  negligibly  small,  and  both 
gas  and  particle  flow  fields  for  one-way  and  two-way 
coupling  are  identical.  Particle  speed  is  practically 
constant  at  1,350  m/s,  and  does  not  change  with  dis¬ 
tance  from  the  nozzle.  Particle  temperature  decreases 
only  due  to  the  radiation  cooling.  When  the  radia¬ 
tion  cooling  is  turned  off,  it  stays  at  about  2,300  K 
throughout  the  plume.  This  is  illustrated  in  Fig.  16, 
where  the  particle  temperature  with  and  without  the 
radiation  cooling  is  shown  in  XY  plane.  Note  here 
that  the  effect  of  this  cooling  on  gas  molecules  has 
not  been  incorporated  in  the  present  model,  and  may 
have  some  impact  on  gas  properties. 

The  distribution  of  particle  number  density  in  XY 
plane  is  presented  in  Fig.  17.  Similar  to  the  flow  inside 
the  nozzle,  at  each  axial  station  downstream  from  the 
nozzle  exit,  the  particle  density  has  a  maximum  off 
the  nozzle  centerline,  which  is  related  to  the  above 
mentioned  fact  that  inside  the  converging  part  of  the 
nozzle  particles  concentrate  in  the  regions  close  to  the 
surface.  Note  also  that  at  about  10  m  downstream 
from  the  nozzle  exit,  the  particle  density  drops  over 
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three  orders  of  magnitude,  proportionally  decreasing 
the  collision  rate  of  gas  molecules  with  particulates. 

This  is  illustrated  in  Fig.  19,  where  the  mean  free 
path  of  gas  molecules  between  collisions  with  particles 
is  shown.  The  mean  free  path  sharply  increases  with 
distance  from  the  nozzle,  and  reaches  about  3  km  at 
9  m  downstream  from  the  nozzle  exit. 

Let  us  now  consider  the  possibility  of  the  interac¬ 
tion  of  high-speed  free  stream  molecules  with  plume 
particulates.  The  number  density  of  the  free  stream 
nitrogen  is  presented  in  Fig.  19.  It  is  clear  that  the 
plume  number  density  in  the  first  9  m  from  the  nozzle 
exit  is  too  high  for  free  stream  molecules  to  penetrate 
to  the  nozzle  centerline.  There  is  a  weak  shock  ob¬ 
served  at  the  wind  side  of  the  plume,  formed  by  the 
free  stream  molecules,  with  the  density  almost  three 
times  as  high  as  in  the  free  stream,  but  there  are  no 
free  stream  molecules  in  the  plume  core  flow  where 
the  particles  reside  (compare  with  Fig.  18).  For  dis¬ 
tances  along  the  nozzle  centerline  larger  than  consid¬ 
ered,  the  collisions  between  the  free  stream  and  plume 
particulates  may  occur,  but  the  particulate  density  is 
expected  to  be  too  small  for  these  collisions  to  have  a 
significant  effect. 

Figures  20  and  21  show  the  total  gas  number  den¬ 
sity  and  translational  temperature,  respectively.  The 
impact  of  particulates  is  visible  in  the  corefiow.  As 
expected,  the  maximum  of  translational  temperature 
is  observed  in  the  plume-atmosphere  interaction  re¬ 
gion.  There  is  also  a  small  local  maximum  at  the 
nozzle  axis  due  to  gas-particle  interaction  inside  the 
nozzle  and  in  a  small  vicinity  downstream  from  the 
nozzle  exit.  The  computations  also  showed  that  with 
the  complete  vibrational  energy  accommodation  the 
vibrational  temperature  of  the  plume  species  is  about 
1,500  K  in  the  plume  coreflow,  and  about  1,000  K 
in  the  particle  free  regions  of  the  plume  and  in  the 
plume-interaction  regions. 

7  Conclusions 

Numerical  investigation  of  a  two-phase  plume  from 
a  small  side  thruster  interacting  with  atmosphere  at 
120  km  was  performed  using  a  combined  continuum- 
kinetic  approach.  The  solution  of  Navier-Stokes  equa¬ 
tions  inside  the  nozzle  was  used  to  specify  a  starting 
surface  for  a  2D  DSMC  computation  of  the  plume  near 
field,  that  was  then  used  for  a  subsequent  3D  DSMC 
computation  of  the  plume-atmosphere  interaction  re¬ 
gion.  The  details  of  the  DSMC  implementation  of  the 
two-way  gas-particulate  coupling  are  presented.  The 
current  implementation  uses  molecular  fluxes  to  cal¬ 
culate  the  number  of  gas-particulate  collisions,  and  is 
based  on  a  statistical  approach  to  calculate  deflection 
angles. 
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The  study  of  sensitivity  of  computational  results 
to  various  parameters  of  the  approach  has  been  per¬ 
formed.  It  was  found  that 

•  For  small  thrusters,  accurate  modeling  of  partic¬ 
ulates  inside  the  converging  part  of  the  nozzle  is 
critical  since  it  practically  determines  the  particle 
divergence  from  the  nozzle  axis  inside  the  plume 

•  The  interaction  between  gas  and  particulates 
(two-way  coupling)  has  an  important  effect  on 
the  flow  in  the  diverging  part  of  the  nozzle 

•  Particles  are  concentrated  in  the  coreflow,  and 
the  divergence  angle  of  particulates  from  the  noz¬ 
zle  centerline  is  about  6  deg 

•  The  kinetic  and  continuum  solutions  of  the  flow 
inside  the  nozzle  are  in  reasonable  agreement  in 
terms  of  gas  as  well  as  particulate  properties 

•  A  40  percent  reduction  in  particle  mass  increases 
the  divergence  angle  in  the  diverging  part  of  the 
nozzle  by  a  factor  of  two 

•  Accommodation  of  vibrational  energy  of  gas 
molecules  at  the  particulate  surface  has  little  ef¬ 
fect  on  vibrational  temperature 

•  The  two-way  coupling  has  negligible  impact  on 
the  flow  after  a  few  centimeters  from  the  nozzle 
exit 

•  Particle  velocities  do  not  change  noticeably  after 
the  nozzle  exit;  particle  temperature  in  the  plume 
decreases  only  due  to  radiative  cooling 

•  There  is  no  interaction  between  the  free  stream 
molecules  and  particulates  in  the  first  ten  meters 
from  the  nozzle  exit;  little  effect  is  expected  from 
this  interaction  in  the  far  field  of  the  plume. 
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Figure  1:  Particle- molecule  collision  in  the  local  coor¬ 
dinate  system. 


Figure  3:  Surface  temperature  distribution  of  parti¬ 
cles  introduced  in  uniform  gas  flow. 
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Figure  4:  Gas  overall  temperature  (K)  field  inside  the 
nozzle  obtained  with  the  continuum  (top)  and  kinetic 
(bottom)  approaches. 
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Figure  5:  Particle  temperature  (K)  field  inside  the 
nozzle  obtained  with  the  continuum  (top)  and  kinetic 
(bottom)  approaches. 
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Figure  7:  Particle  density  (kg/m3)  field  inside  the 
nozzle  obtained  with  the  continuum  (top)  and  kinetic 
(bottom)  approaches. 
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Figure  6:  Particle  axial  velocity  (m/s)  field  inside  the 
nozzle  obtained  with  the  continuum  (top)  and  kinetic 
(bottom)  approaches. 


Figure  8:  Impact  of  two-way  coupling  on  gas  axial  ve¬ 
locity  (m/s)  inside  the  nozzle.  Top,  one-way  coupling; 
bottom,  two-way  coupling. 
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Figure  9:  Impact  of  two-way  coupling  on  particle  axial 
velocity  (m/s)  inside  the  nozzle.  Top,  one-way  cou¬ 
pling;  bottom,  two-way  coupling. 


Figure  11:  Impact  of  particle  mass  on  gas  axial  veloc¬ 
ity  (m/s)  inside  the  nozzle.  Top,  aluminum  particles; 
bottom,  alumina  particles. 
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Figure  10:  Impact  of  particle  mass  on  particle  axial 
velocity  (m/s)  inside  the  nozzle.  Top,  aluminum  par¬ 
ticles;  bottom,  alumina  particles. 


Figure  12:  CO  vibrational  temperature  (K)  inside  the 
nozzle  for  different  vibrational  energy  accommodation 
coefficients  on  particle  surface.  Top,  full  accommoda¬ 
tion;  bottom,  no  accommodation 
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Figure  13:  Impact  of  two-way  coupling  on  gas  axial 
velocity  (m/s)  in  the  plume  near  field.  Top,  one-way 
coupling;  bottom,  two-way  coupling. 


Figure  14:  Impact  of  two-way  coupling  on  gas  transla¬ 
tional  temperature  (K)  in  the  plume  near  field.  Top, 
one-way  coupling;  bottom,  two-way  coupling. 
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Figure  15:  Particle  temperature  (K)  and  streamlines 
in  the  plume  near  field. 
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Figure  16:  Particle  temperature  (K)  in  the  plume- 
atmosphere  interaction  region  with  radiation  cooling 
turned  off  (top)  and  on  (bottom).  Here,  white  region 
shows  the  location  of  the  starting  surface  for  3D  mod¬ 
eling 
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Figure  17:  Particle  number  density  (1/m3)  in  the 
plume-atmosphere  interaction  region. 
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Figure  18:  Gas  mean  free  path  between  collisions  with  Figure  20:  Gas  number  density  (1/m3)  in  the  plume- 

particles  (m)  in  the  plume-atmosphere  interaction  re-  atmosphere  interaction  region 

gion. 


6. 439**16 
l*288e-t-17 
1  *  832*1-17 
2+576e*17 
3 *220**17 
3*864e*17 
4*507**17 
5,151**17 
5,795**17 
G, 438**1 7 

7*083ei 17 

7*727e+17 
8*37le*17 
8*0l5e*l7 
3.658**17 
1  *  O3Oc+10 
1*036**16 
1,158**16 
1*223**18 
1*288**18 
1* 352**1 8 
1*417**18 
1*481**18 
1,545**18 


Figure  19:  Free  stream  nitrogen  number  density  Figure  21:  Gas  translational  temperature  (K)  in  the 
(1/m3)  in  the  plume-atmosphere  interaction  region.  plume- atmosphere  interaction  region. 


